Se utilizará los datos de censo para 8 condados centrales del estado de New York
ny <- st_read('data/NY8_utm18.shp')
## Reading layer `NY8_utm18' from data source
## `C:\Users\Usuario\Documents\Universidad\2022-I\Espacial\ArealData\data\NY8_utm18.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## Projected CRS: WGS 84 / UTM zone 18N
plot(st_geometry(ny), axes = T)
Un archivo gal es una matriz de pesos o vecinos
producida por el software GeoDa.
library(spdep)
NY_nb <- read.gal('data/NY_nb.gal', # read.gal carga el archivo de clase nb
region.id = as.numeric(row.names(ny))-1) # Debe empezar en 0
Syracuse <- filter(ny, AREANAME == 'Syracuse city')
Sy0_nb <- subset(NY_nb, ny$AREANAME == 'Syracuse city')
Se realiza un subset del área de trabajo por temas de simplicidad.
plot(st_geometry(Syracuse), border = 'grey60', axes = T)
plot(Sy0_nb, st_geometry(Syracuse), pch = 19, cex = 0.6, add = T)
Las relaciones entre n observaciones vecinas (regiones) son representadas por un objeto de clase nb. Estos objetos son una lista de longitud n con el índice de la región y un vector de enteros con los vecinos de cada región (índices). Si una región (polígono) no tiene vecinos, el componente tendrá un 0.
El objeto de clase nb (NY_nb) tiene incluidos métodos de
print, summary y plot entre
otros.
Summary presenta una tabla de la distribución del número de enlaces y reporta asimetría y presencia de observaciones no vecinas.
summary(Sy0_nb)
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 346
## Percentage nonzero weights: 8.717561
## Average number of links: 5.492063
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9
## 1 1 5 9 14 17 9 6 1
## 1 least connected region:
## 164 with 1 link
## 1 most connected region:
## 136 with 9 links
La asimetria está presente cuando i es un vecino de j, pero j no es un vecino de i.
El resultados de la función summary brinda información
acerca de:
Se observa también la distribución de las relaciones que indican que el número máximo de vecinos que se encontraron fueron 9 y estuvieron con respecto a un poligono, 17 de las 63 regiones presentes tienen 6 uniones, etc.
card(Sy0_nb)
## [1] 5 5 2 7 7 5 4 5 7 5 5 7 5 7 6 6 6 6 6 4 4 7 6 8 4 3 3 9 6 8 6 6 6 8 6 4 4 8
## [39] 6 6 7 5 7 6 4 5 3 5 6 4 6 7 4 8 5 1 5 8 5 5 6 3 3
La función card permite obtener de forma rápida el
número de vecinos de cada región.
Se crearán objetos de lista de vecinos para uso posterior en el
documento, para esto se utiliza la función knearneigh para
obtener el índice de los puntos de k vecinos cercanos, esta
función recibe como parametro de ingreso la lista de coordenadas de los
puntos y un valor k que indica el número de puntos vecinos
cercanos que se desean retornar. El resultado de la operación anterior
se ingresa a la función knn2nb para convertir el objeto
resultante en una lista de vecinos de la clase nb.
coords <- st_point_on_surface(st_geometry(Syracuse))
ids <- row.names(Syracuse)
Sy8_nb <- knn2nb(knearneigh(coords, k = 1), row.names = ids)
Sy8_nb
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 63
## Percentage nonzero weights: 1.587302
## Average number of links: 1
## Non-symmetric neighbours list
la función nbdists obtiene las distancias entre los
puntos de los objectos vecinos calculados en el paso anterior. Esto es
util porque con ella se obtiene la distancia minima necesaria
para asegurar que todas las áreas tengan al menos un
vecino.
Posteriormente se utiliza la función dnearneigh para
obtener la lista de objetos vecinos (nb) que se encuentran a una
determinada distancia (d1, d2).
dsts <- unlist(nbdists(Sy8_nb, coords))
Sy11_nb <- dnearneigh(x = coords, d1 = 0, d2 = 0.75*max(dsts))
Sy11_nb
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 156
## Percentage nonzero weights: 3.930461
## Average number of links: 2.47619
## 8 regions with no links:
## 1 10 20 36 46 57 60 61
Como resultado, se observa que existen 8 regiones que no poseen vecinos (1, 10, 20…).
Los pesos espaciales se pueden ver como una lista de pesos indexados por una lista de vecinos, donde el peso de la relación entre i e j es el k-ésimo elemento del i-ésimo componente de la lista de pesos.
La función nb2listw toma un objeto de lista de vecinos
(nb) y lo convierte en un objeto de pesos. La conversión por defecto es
W, donde los pesos varian entre la unidad dividido por el más
largo y el más pequeño número de vecinos y la suma de los pesos para
cada entidad de área es la unidad.
Sy0_lw_w <- nb2listw(neighbours = Sy0_nb)
Sy0_lw_w
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 346
## Percentage nonzero weights: 8.717561
## Average number of links: 5.492063
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 63 3969 63 24.78291 258.564
Otro tipo de conversión se obtiene al configurar style = B, con lo cual se obtiene un peso de 1 por cada relación vecina, pero en este caso, la suma de los pesos para áreas difieren de acuerdo al número de vecinos que las áreas tengan.
Se recomienda utilizar pesos binarios cuando se tiene poco conocimiento acerca del proceso que se está estudiando.
try(
expr = {
nb2listw(neighbours = Sy11_nb, style = 'B', zero.policy = F)
}
)
## Error in nb2listw(neighbours = Sy11_nb, style = "B", zero.policy = F) :
## Empty neighbour sets found
Cuando existen regiones sin vecinos es necesario especificar el
argumento zero.policy = T para que no arroje error.
Sy0_lw_D1 <- nb2listw(neighbours = Sy11_nb, style = 'B', zero.policy = T)
print(Sy0_lw_D1, zero.policy = T)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 156
## Percentage nonzero weights: 3.930461
## Average number of links: 2.47619
## 8 regions with no links:
## 1 10 20 36 46 57 60 61
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 55 3025 156 312 2240
Hay dos tipos de test de autocorrelación espacial que se aplican a este tipo de datos:
Es la forma más común de abordar la autocorrelación espacial. Se calcula como la razón del producto de la variable de interés y su rezago espacial, con el producto cruzado de la variable de interés y ajustada por los pesos espaciales usados.
\[I = \frac{n}{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}}\ \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} (y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^{n} (y_i - \bar{y})^2}\] donde:
Centrarse en la media es equivale a afirmar que el modelo correcto tiene una media constante, y que cualquier patrón restante después del centrado es causado por las relaciones espaciales codificadas en los pesos espaciales.
Para esto, se calcula la esperanza como:
\[E(I) = \frac{-1}{n-1}\]
El siguiente paso es tomar el valor observado y restarlo al valor esperado analíticamente (I - E(I)) para finalmente dividirlo por la raíz cuadrada de la varianza analítica. El resultado es comparado con la distribución normal para encontrar el valor de probabilidad del estadístico observado bajo la hipótesis nula de ausencia de dependencia espacial para los pesos espaciales elegidos. Los resultados pueden depender de la elección de las ponderaciones, de los vecindarios y de que los supuestos sean cumplidos.
\[H_0: Ausencia \ de \ dependencia \ espacial\] \[H_1: El \ estadístico \ observado \ es \ significativamente \ más \ grande \ que \ el \ valor \ esperado\]
El dominio del I de Moran es [-1, 1]. Se evidencia autocorrelación espacial positiva cuando I > 0. Esto quiere decir que las unidades de análisis vecinas tienden a ser similares.
moran.test(x = ny$Cases, listw = nb2listw(NY_nb))
##
## Moran I test under randomisation
##
## data: ny$Cases
## weights: nb2listw(NY_nb)
##
## Moran I statistic standard deviate = 3.9778, p-value = 3.477e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.146882990 -0.003571429 0.001430595
Por defecto, el moran.test usa el supuesto de
aleatorización lo que difiere del supuesto de normalidad más simple al
introducir un término de corrección basado en la curtosis de la variable
de interés. Cuando la curtosis corresponde a la de una variable con
distribución normal, las dos suposiciones arrojan la misma varianza,
pero a medida que la variable se aleja de la normalidad, la suposición
de aleatorización compensa aumentando la varianza y disminuyendo la
desviación estándar.
moran.test(x = ny$Cases, listw = nb2listw(NY_nb), randomisation = F)
##
## Moran I test under normality
##
## data: ny$Cases
## weights: nb2listw(NY_nb)
##
## Moran I statistic standard deviate = 3.9731, p-value = 3.547e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.146882990 -0.003571429 0.001433975
En este caso, se observa que la varianza no cambia, por lo tanto, el supuesto de normalidad se cumple.
Es posible también realizar test de Monte Carlo para evaluar la hipótesis de dependencia espacial con el test de Moran I.
moran.mc(x = ny$Cases, listw = nb2listw(NY_nb), nsim = 999)
##
## Monte-Carlo simulation of Moran I
##
## data: ny$Cases
## weights: nb2listw(NY_nb)
## number of simulations + 1: 1000
##
## statistic = 0.14688, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
Conociendo las medidas globales de asociación espacial, es posible descomponer estas medidas con el fin de construir test locales que buscan detectar clusters (Observaciones con vecinos muy similares) y hotspots (Observaciones con vecinos muy diferentes), además se pueden utilizar para conocer el impacto de los estadísticos individuales en sus contrapartes globales.
El gráfico de dispersión de Moran por convención posiciona la variable de interés en el eje x y los valores de rezago espacial (La suma de valores de los vecinos espacialmente ponderada) en el eje y. El índice global de Moran es una relación lineal entre estos valores y es dibujado como la pendiente.
(ggplot(moran_plot, aes(x, wx))+
geom_point(aes(shape = is_inf, color = is_inf), size = 0.75) +
scale_color_manual(values = c("black", "red")) +
geom_hline(yintercept = moran_plot$wx %>% mean(),
linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = moran_plot$x %>% mean(),
linetype = "dashed", alpha = 0.5) +
geom_text(aes(x, wx, label = labels),
data = filter(moran_plot, is_inf),
size = 2.5) +
geom_smooth(method = "lm", se = F, size = .6, col = "black") +
labs(x = "Cases", y = "Spatial lagged cases",
title = "Moran Scatterplot",
shape = "Significativo", color = "Significativo") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))) %>%
plotly::ggplotly()
El gráfico se divide en 4 cuadrantes. El cuadrante de la esquina superior derecha indica áreas que tienen un alta incidencia de casos de leucemia y están rodeadas por otras áreas que tienen nivel superior al promedio de casos (high-high). En la esquina inferior izquierda se encuentran áreas con baja incidencia de casos de leucemia rodeadas por áreas con incidencia de leucemia por debajo del promedio (low-low). Las áreas high-high son conocidas como hotspots y las áreas low-low como coldspots. En la diagonal opuesta se tienen valores atípicos espaciales, que indica que son áreas que están rodeadas por vecinos muy diferentes a ellos.
Los valores locales del índice de Moran se construyen como los n componentes sumados para alcanzar el índice de Moran global.
\[I_i = \frac{(y_i - \bar{y}) \sum_{j=1}^n w_{ij}(y_j - \bar{y})}{\frac{\sum_{i = 1}^n (y_i - \bar{y})^2}{n}}\]
Al igual que en el caso de la estadística global, se puede comprobar la divergencia de los estadísticos locales de los valores esperados, bajo los supuestos de normalidad y aleatoriedad analíticamente y usando métodos como el de punto de silla y exacto.
local <- as.data.frame(localmoran(ny$Cases, listw = nb2listw(NY_nb, style = "C")))
local2 <- as.data.frame(localmoran.sad(lm(Cases ~ 1, ny), nb = NY_nb, style = "C"))
local3 <- as.data.frame(localmoran.exact(lm(Cases ~ 1, ny), nb = NY_nb, style = "C"))
ny%>%
select(Cases)%>%
mutate(p_value = moran_plot$is_inf,
rezago = lag.listw(nb2listw(NY_nb, style = "B"), Cases),
cuadrante = case_when(
(Cases > mean(Cases)) & (rezago > mean(rezago)) &
(p_value == T) ~ 'High-High',
(Cases <= mean(Cases)) & (rezago <= mean(rezago)) &
(p_value == T) ~ 'Low-Low',
(Cases > mean(Cases)) & (rezago <= mean(rezago)) &
(p_value == T) ~ 'High-Low',
(Cases <= mean(Cases)) & (rezago > mean(rezago)) &
(p_value == T) ~ 'Low-High',
TRUE ~'NS'
),
index = moran_plot$labels) -> resultados
resultados_coords <- filter(resultados, p_value == T) %>%
mutate(x = (st_centroid(resultados$geometry) %>% st_coordinates())[1],
y = (st_centroid(resultados$geometry) %>% st_coordinates())[2])
(ggplot(resultados)+
geom_sf(aes(fill = cuadrante), alpha = 0.8)+
scale_fill_manual(values = viridisLite::viridis(5))+
labs(title = 'Lugares con influencia de Leucemia \n',
x = 'Longitud', y = 'Latitud')+
geom_sf_text(data = resultados_coords,
mapping = aes(x, y, label = index),
color = "white") +
theme_bw() +
theme(plot.title = element_text(size = 14, face = 'bold', hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank())) %>%
plotly::ggplotly()
Se puede calcular una evaluación de riesgo constante asumiendo una distribución Poisson con un valor de \(\lambda\) igual a la probabilidad de ocurrencia de la enfermedad, se estandariza la variable tomando la propiedad de que tanto la media como la desviación estandar son iguales en esta distribución, se calculan los rezagos y se compara el indice de morán actual con el de riesgo constante.
ny %>%
transmute(p_afectada = sum(Cases)/sum(POP8)*POP8,
casos_std = (Cases - p_afectada)/sqrt(p_afectada),
rezago = stats::lag(nb2listw(NY_nb, style = 'C'), casos_std),
indice_riesgo = casos_std * rezago,
indice = local$Ii) -> riesgo
m_riesgo <- ggplot(riesgo)+
geom_sf(aes(fill = indice_riesgo))+
labs(title = 'Mapa de riesgo Constante')+
scale_fill_gradientn(colours = viridisLite::viridis(20, begin = 0.17,
end = 0.98),
limits = c(-2.2, 3.7)) +
theme_bw()+
theme(plot.title = element_text(size = 14, face = 'bold', hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank())
m_actual <- ggplot(riesgo)+
geom_sf(aes(fill = indice))+
labs(title = 'Mapa de riesgo Standard')+
scale_fill_gradientn(colours = viridisLite::viridis(20, begin = 0.17,
end = 0.98),
limits = c(-2.2, 3.7)) +
theme_bw() +
theme(plot.title = element_text(size = 14, face = 'bold', hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank())
cowplot::plot_grid(m_actual, m_riesgo, align = 'hv', ncol = 2)
# set.seed(1234)
# nsim <- 999
# N <- length(rni)
# sims <- matrix(0, ncol = nsim, nrow = N)
# for (i in 1:nsim) {
# y <- rpois(N, lambda = rni)
# sdCRi <- (y - rni)/sqrt(rni)
# wsdCRi <- lag(lw, sdCRi)
# sims[, i] <- sdCRi * wsdCRi
# }
# xrank <- apply(cbind(I_CR, sims), 1, function(x) rank(x)[1])
# diff <- nsim - xrank
# diff <- ifelse(diff > 0, diff, 0)
# pval <- punif((diff + 1)/(nsim + 1))
monte_carlo_sim <- function(nsim, prop, Ii, nb){
N <- length(prop)
sims <- matrix(0, ncol = nsim, nrow = N)
for (i in 1:nsim) {
y <- rpois(N, lambda = prop)
z <- (y - prop)/sqrt(prop)
rezago <- stats::lag(nb2listw(nb, style = "C"), z)
sims[, i] <- z * rezago
}
xrank <- apply(cbind(Ii, sims), 1, function(x) rank(x)[1])
diff <- nsim - xrank
diff <- ifelse(diff > 0, diff, 0)
pval <- punif((diff + 1)/(nsim + 1))
return(pval)
}
mc_sim <- riesgo %>%
select(p_afectada, indice_riesgo) %>%
mutate(pvalue = monte_carlo_sim(999, p_afectada,
indice_riesgo, NY_nb))
ggplot(mc_sim)+
geom_sf() +
geom_sf(aes(fill = pvalue),
data = filter(mc_sim, pvalue < 0.1 | pvalue > 0.9))+
labs(title = 'Mapa de riesgo aleatorio')+
scale_fill_gradientn(colours = viridisLite::viridis(20, begin = 0.17,
end = 0.98)) +
theme_bw() +
theme(plot.title = element_text(size = 14, face = 'bold', hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank())